clear
clc

inf=zeros(68,2);

parfor zz=1:68
tic    
%%%%%loading data    
w=load('competitions.mat');
COMP=w.comp_list(zz,1);
x=load('data.mat');
y=load('qestimates.mat');
competitionid=x.competitionid;
maxscore = x.maxscore(competitionid==COMP,1);
pubscore_normal = x.pubscore_normal(competitionid==COMP,1);
length_days= x.length_days(competitionid==COMP,1);
length_days=length_days(1);
rewardquantity= x.rewardquantity(competitionid==COMP,1);
rewardquantity=rewardquantity(1);
mergers_data= x.totmergers(competitionid==COMP,1);
mergers_data=mergers_data(1);
nSim = 3000;

%%%%%state space parameters
N=40; %number of individuals
T=360; %number of periods
Svalues=unique([unique(maxscore);max(maxscore)+[0.002:0.002:0.08]']);
nScore = 20;
Svalues = Svalues(3:nScore+2,1);
nsubs_data=size(maxscore,1);
lambda1=0.75;
lambda2=1-lambda1;

%%%%%q-function
betachangehat=[y.beta0(y.cid==COMP,1) ; y.beta1(y.cid==COMP,1); y.beta1tr(y.cid==COMP,1)];
probchange = @(score) (exp(betachangehat(1) + betachangehat(2)*score)./(1+exp(betachangehat(1) + betachangehat(2)*score)));
q_sp=probchange(Svalues);
%q_sp(end,1)=0;
probchange = @(score) (exp(betachangehat(1) + betachangehat(3) + betachangehat(2)*score)./(1+exp(betachangehat(1)+ betachangehat(3) + betachangehat(2)*score)));
q_team=probchange(Svalues);
%q_team(end,1)=0;

%load estimates
asdf=load(sprintf('Estimates/%s_%02d.mat', 'estimates',COMP));
sigmaS=asdf.sigmaS;
sigmaT=asdf.sigmaT;

%%%%%standard errors
[mergers,subs,maxs,Omega]=equilibrium(sigmaS,sigmaT,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
G=zeros(2,2);
epsilonT=1e-2;
epsilonS=1e-4;
L=[(mergers-mergers_data)/mergers_data;(subs-nsubs_data)/nsubs_data];
[m1,s1,~,~]=equilibrium(sigmaS,sigmaT+epsilonT,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
L1=[(m1-mergers_data)/mergers_data;(s1-nsubs_data)/nsubs_data];
G(:,1)=[L1-L]/(epsilonT);
[m1,s1,~,~]=equilibrium(sigmaS+epsilonS,sigmaT,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
L1=[(m1-mergers_data)/mergers_data;(s1-nsubs_data)/nsubs_data];
G(:,2)=[L1-L]/(epsilonS);
VarCov=(inv(G'*G)*G'*(L*L')*G*inv(G'*G));
StErrorSigma=diag(VarCov/nsubs_data).^0.5;
inf(zz,:) = StErrorSigma';
end

clearvars -except inf
save sigma_inference

